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Abstract 

We analyze a piecewise-linear FitzHugh-Nagumo model. The system exhibits 
a canard near which both small amplitude and large amplitude periodic orbits 
exist. The addition of small noise induces mixed- mode oscillations (MMOs) in 
the vicinity of the canard point. We determine the effect of each model parameter 
on the stochastically driven MMOs. In particular we show that any parameter 
variation (such as a modification of the piecewise-linear function in the model) that 
leaves the ratio of noise amplitude to time-scale separation unchanged typically has 
little effect on the width of the interval of the primary bifurcation parameter over 
which MMOs occur. In that sense, the MMOs are robust. Furthermore we show 
that the piecewise-linear model exhibits MMOs more readily than the classical 
FitzHugh-Nagumo model for which a cubic polynomial is the only nonlinearity. 
By studying a piecewise-linear model we are able to explain results using analytical 
expressions and compare these with numerical investigations. 



1 Introduction 

Oscillatory dynamics involving oscillations with greatly differing amplitudes, known as 
mixed- mode oscillations (MMOs), see Fig. [H are important in neuron models [1] and in 
a multitude of chemical reactions (2j |3], refer to |4] for a recent review. Yet there are 
many open questions regarding the creation, robustness and bifurcations of MMOs. A 
variety of mechanisms generate MMOs in deterministic systems. Alternatively MMOs 
may be noise-induced; there are also several scenarios by which this may occur. 

*The authors acknowledge support from an NSERC Discovery Grant. 
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We study the following form of the FitzHugh-Nagumo (FHN) model with small, 
additive, white noise: 

dv = (f(v) -w)dt, 

dw = e(av - aw - A) dt + D dW , 

where v represents a potential, w is a recovery variable and If is a standard Brownian 
motion. The FHN model is used as a prototypical model of excitable dynamics in a range 
of scientific fields [5j [6]. Here a is a positive constant and A £ R, which is regarded as 
the main bifurcation parameter, controls the growth of oscillations, as seen below. The 
small parameter £ < 1, represents the time-scale separation and D <C 1 is the noise 
amplitude (e, D > 0). Values of e and D used in, for instance |8|, are no larger than 
the values considered here. By scaling we may assume a = 1, except in the special case 
a = which corresponds to the van der Pol model (and in this case we may further 
assume a = 1). We assume that / : R — > R is continuous and roughly of cubic shape. 
For simplicity, we assume that / has a local minimum at (0, 0) and a local maximum at 
(1,1), regardless of the precise function chosen. 

If / is a cubic, as originally taken by FitzHugh [9] and Nagumo et. al. [10], then, by 
the above requirements, the cubic must be 

f( v ) = 3v 2 - 2v 3 . (2) 

Fig. [2]-A illustrates the role of the parameter A for (OQ) with ([2]) in the absence of noise. 
A small amplitude periodic orbit is created in a Hopf bifurcation at A = 0. For the 
parameters used in Fig. [2], this periodic orbit is stable and its amplitude increases with A. 
Near A c the amplitude increases exponentially. This rapid growth is known as a canard 
explosion and is due to time-scale separation and global dynamics [HI [T5J [131 EH] ■ The 
value of the canard point, A c , which is well-defined for smooth systems [JHHTB], decreases 
to zero with e, as shown in Fig. [3]-A. Over an order e range of A values, ([1]) with (J2]) 
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Figure 1: A time series illustrating MMOs exhibited by ([T]) with ([3]). The parameter 
values are the same as in Fig. |4]A. 



2 



smooth 



1.5 



1- 



0.5 



1.5 



B 



PWL 



0.5 



0.02 



0.04 




Figure 2: Bifurcation diagrams of ([T]) in the absence of noise (i.e. D = 0) with (J2J) in 
panel A and with ([3]) in panel B. In each panel the solid curve for A > corresponds to 
the maximum v- value of a stable periodic orbit; the remaining curves correspond to the 
equilibrium which is unstable for A > 0. In panel A a canard explosion occurs near the 
canard point, A c ; in panel B a canard explosion occurs near Ai at which point the stable 
periodic orbit has a maximum value of 1. The parameter values used are e = 0.04, 
(a,a) = (4,1), (r)L,VR) = ("2,-1) and (v hWl ) = (0.1,0.05). 



may either settle to equilibrium, exhibit small amplitude oscillations, or exhibit large 
amplitude oscillations (relaxation oscillations). 

As in [TTt IT8] . here we study a piecewise-linear (PWL) FHN model so that, in the 
presence of noise, the system is amenable to a rigorous analysis without the need for an 
approximation or limiting scenario. PWL models are commonly used in circuit systems 
[I~9 t 120 1 |2T] . A PWL version of a driven van der Pol oscillator is studied in [22] to explain 
the breakdown of canards in experiments. We consider the continuous, PWL function 

TJlV , v < 

f(r) = { mv r ^ °V-J^ > ( 3 ) 

' r] 2 {v — v-i) + Wx , Vx < v < 1 w 
rj R (v - 1) + 1 , v > 1 

where < V\, Wi < 1, t]l, Vr < 0, and 

wi 1-wx 
Vi = — , V2 = • (4) 

Vt 1-Vx 

We state the particular form here in order to briefly illustrate key differences between the 
smooth and PWL FHN models. Further motivation for ([3]) is given in §|2] and shown in 
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Figure 3: Two parameter bifurcation diagrams of the smooth and PWL versions of ([1]) 
with the same parameter values as in Fig. |2j The smooth system has a well-defined 
canard point, A c [151 US]) whereas for the PWL system we consider the two values, X V1 
and Ai, described in the text. In both panels we have indicated the attracting solution 
for each region bounded by the solid curves. The dotted curve in panel B corresponds to 
the approximation ( 1T4"j) derived below; £ crit is given by (ITTi) . Note that in contrast to the 
remainder of this paper, in panel A the distinction between small and large oscillations 
is determined by A c and not fl5]). 



Fig. [51 As shown in Fig. |2]T3, ([T}) with ([3]) may exhibit a canard explosion. The canard 
point, A c , is not well-defined for this system because it lacks global differentiability. 
Instead we consider the values, X Vl and A 1; at which the maximum v- value of the periodic 
orbit of ([1}) with ([3]) in the absence of noise is v\ and 1 respectively. The piecewise 
nature of fl3]) leads to a natural classification of periodic orbits and oscillations of (JTJ) 
with ([3]). (Typically we refer to one complete revolution about the equilbrium as a single 
oscillation.) With Fig. |2]T3 in mind, if t> max is the maximum v- value of a periodic orbit 
or single oscillation we declare that the orbit or oscillation is 

small if < i> max < v i , 
medium if v± < i> max < 1 , (5) 
large if v m3X > 1 . 

Fig. [3]T3 illustrates typical dependence of X Vl and Ai on e. In particular we notice that 
for a fixed choice of the slopes, rjj, in (E]), the PWL version of the FHN model does 
not exhibit small oscillations for arbitrarily small e. This is because the two eigenvalues 
associated with the equilibrium for small A > are real-valued for sufficiently small e 
negating the possibility of small oscillations, see §|2j Unlike for the van der Pol model, 
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Figure 4: A trajectory of ([1]) with ([2]) in panel A, and ([1]) with ([3]) in panel B. The 
parameter values used are the same as in Fig. |2j with also A = 0.028, D = 0.0008. For 
these parameter values both the smooth and PWL models are tuned to near the canard 
explosion. 



values of e that are relevant for the FHN model are usually sufficiently large for small 
oscillations to be important in the PWL model. 

The effect of noise in (TjQ) has seen significant recent attention, see for instance [23 | l2"3 j 
[25] . Noise may induce regular oscillations in (CQ) when in the absence of noise there are 
no oscillations. There is more than one mechanism that may cause this, most notably 
stochastic resonance [26] (when a small periodic forcing term is present in addition to 
noise), coherence resonance [25] (usually when the system is quiescent in the absence of 
noise), and self -induced stochastic resonance |27J (involving relatively large noise that 
drives oscillations of periods different from that of the deterministic system). 

If A is tuned to values near the canard explosion, in the presence of noise the system 
may exhibit both small amplitude and large amplitude oscillations, i.e. MMOs, as shown 
in Figs. H] and [TJ Similar MMOs are described in [28] for very small noise by a careful 
choice of parameter values. For a version of ([1]) that contains nonlinearity in the w 
equation to better mimic neural behaviour, it has been observed that when A is chosen 
to be just prior to the canard point the frequency of relaxation oscillations increases with 
noise amplitude [8]. Noise- induced MMOs have been described for three coupled FHN 
systems near a canard [29]. A signal-to-noise ratio may be defined to quantitatively 
determine dominant frequencies [30] . Noise- induced MMOs may arise via a different 
mechanism in the case that the Hopf bifurcation is subcritical [31]. There are a variety 
of mechanisms for MMOs in three-dimensional systems that we do not consider here, 
see for instance [32J ES] arid references in [1] . 
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In this paper we study noise-driven MMOs in ([T]) with (j3j. We use analytical methods 
to identify parameter values for which MMOs occur and describe the dependence of each 
model parameter on MMOs. For typical values of the noise amplitude, D, MMOs occur 
over some interval of positive A-values. In order to find such intervals we determine 
exit distributions for forward orbits of §Q with through various cross-sections of 
phase space. The exit distributions allow us to deduce the amplitude of oscillations and 
consequently find intervals of MMOs. We show that MMOs are robust in the sense that 
large variations in other model parameters can have minimal effect on the width of the 
A-intervals. 

We note that the model we consider has additive noise in the tw-equation only, 
as in, for instance, [8j [25] . This choice allows some simplifications in demonstrating 
the analytical method, while still capturing qualitatively the behavior that would be 
observed for more general additive noise. Throughout the paper we indicate where this 
assumption allows some simplification in the analysis, and we indicate the differences 
that would need to be addressed for the case of noise also in the ^-equation. 

The remainder of the paper is organized as follows. Section [2] briefly overviews PWL 
FHN models and provides an analysis of ([T|) with ([3]) in the absence of noise. Here we 
explain the Hopf-like bifurcation at A = that creates stable oscillations and describe 
equations for X V1 and Ai, Fig. [3J Calculations of exit distributions are detailed in §|3J 
Here we also describe the method by which we use these distributions to find parameter 
values corresponding to MMOs. Section H] combines the analysis of the previous sections 
to determine the effect of each model parameter on MMOs. Finally conclusions are 
presented in §|5j 

2 Properties of the deterministic system 

Analytical results may be derived for (pQ) when f(v) is a PWL function. Arguably 
the simplest continuous, PWL function that one can use for f(v) consists of three line 
segments (one of them being the straight connection between (0,0) and (1,1)). The 
FHN model with this function is well-studied [SJ ES], refer to [3S] for the van der Pol 
system. However, with this three-piece PWL function, (Cp) does not exhibit a canard, as 
shown in [37], and so we do not consider it further. Consequently, as in [221 EB], we use 
two line segments between (0,0) and (1, 1) denoting the intermediate point by (vi,w\) 
and the slopes by rjj, specifically (j3J), as shown in Fig. EJ If instead f{v) contains multiple 
line segments left of (0, 0) such that the slopes of the two lines meeting at (0, 0) are ±77!, 
multiple coexisting attractors commonly exist for small A which leads to complications 
that we do not study here. For simplicity we do not consider f(v) comprised of more 
than four line segments. For canards in PWL FHN models with many segments we refer 
to reader to the recent work of Rotstein et. al. [55] . 

In the absence of noise (i.e. when D = 0), (0Q) with (J3J) is a continuous, two- 
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Figure 5: The nullclines of (JTJ with d3J) for small A > 0. Potential equilibria, (v*,w*), 
lie at the intersection of the nullclines. If rjj < - for every j, then the system has a 
unique equilibrium for all values of A. 



The phase space, 



dimensional, PWL, ordinary differential equation system: 

v = f( v ) ~ w , 

w = e(av — aw — A) . 

is divided into four regions 

TZ L = {(v,w) \ v <0,w ER} , 
TZi = {(v,w) | < v < Vi,w e 1 
TZ 2 — {(v,w) | Vi < v < 1, w e 1 
^ R = {(v,u;) | v > 1,10 e E} , 

0, u = Vi and f 



(6) 



(7) 



1, on which the system is 



by the three switching manifolds, v 
non-differentiable. 

Each linear component of with (J3J) has a unique equilibrium, (v*,w*), see Fig. [5] 
(unless a = ar]j in which case the relevant v- and w-nullclines are parallel). In the 
terminology of piecewise-smooth dynamical systems, each (v*,w*) is either admissible 
(lies in the closure of TZj) or virtual (lies outside the closure of IZj 
and the eigenvalues, pj, associated with each (v*,w*) are 



The Jacobian, Aj, 



A; 



Vj -i 
ea —ea 



P.i 



| (rjj - ea ± a/ (rjj + ea) 2 - Asc^j . 



(8) 
(9) 



We assume 



i] L < —ea — 2y/ea 



a 

ecr < rji < — 



(10) 



such that (v1,w* L ) is an attracting node and (v*,w*) is either a repelling node or a 
repelling focus as determined by the sign of (771 + ea) 2 — Asa. The restriction ffTOj) 
ensures that stable oscillations are created at A = 0, as shown below. 

The bifurcation at A = that results from the interaction of an equilibrium with 
the switching manifold, v = 0, is an example of a discontinuous bifurcation [39., HDJ HE]- 
Effectively, eigenvalues that determine the stability of the admissible equilibrium change 
discontinuously as the equilibrium crosses the switching manifold at A = 0. In general, 
a bifurcation is expected to occur if one or more eigenvalues "jump" across the imagi- 
nary axis at the crossing. Such a bifurcation may be analogous to a smooth bifurcation 
or it may be unique to piecewise-smooth systems [3D]. For two-dimensional systems, 
codimension-one, discontinuous bifurcations involving a single smooth switching mani- 
fold have been completely classified [311 H2] • 

For the PWL system ([6]) with fl3]), an attracting periodic orbit is born at the dis- 
continuous bifurcation, A = 0. The relative size of the periodic orbit for small A > is 
dependent upon whether the equilibrium, (v^, 10*), is a node or a focus. If {v* L , w* L ) is an 
attracting node and (vl,w*) is a repelling node, invariant lines corresponding to eigen- 
vectors prevent the creation of a local periodic orbit corresponding to a small oscillation 
[IT] . The periodic orbit generated at A = has large amplitude (corresponding to a 
relaxation oscillation). Specifically, as A — > + , the maximum value of v of the periodic 
orbit limits on a value greater than 1. 

If instead {v{, w^) is a repelling focus, then the bifurcation is a discontinuous analogue 
of a Hopf bifurcation in that a periodic orbit is created locally. Unlike for a classical 
Hopf bifurcation, the periodic orbit grows in size linearly with respect to A (see Fig. |2]-B) 
which is typical for piecewise-smooth systems. 

The value of e for which the square-root term in (jDJ vanishes is the critical value of e 
(see Fig. [3]-B) above which the periodic orbit created at A = is small and below which 
this orbit is large, and is given by 



The curves A = X Vl (e) and A = Ai(e), Fig. |3]-B, which bound the region of medium os- 
cillations, emanate from (A,e) = (0,£ cr i t ). Since the underlying system is PWL, we may 
obtain analytical expressions relating to these curves by deriving the explicit solution to 
the flow of each linear component of ([I]) with ([3]). We let (v^\t\ v , w ), w^'(t; v , w )) 
denote the solution to the linear component of ODJ with (E]) corresponding to TZj, for an 




(11) 
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arbitrary initial condition, (vq,Wq). For instance: 



v w (t; ^0,^0) 
w w (t; v ,w ) 



cos(wit) + sin^t) 



sin(oJit) 

-^sin(a;it) 
cos(wit) - 2^ sin(o;it) 





Uo-«i 


+ 






_ w - w\ _ 







(12) 



which equals the solution to (EJ) with (|3]) for the same initial condition whenever (i/ 1 ) (t) , itr 1 ) (t)) 
lies in the closure of 1Z\ at all times between and t, and where 



(13) 



Unfortunately we cannot in general explicitly solve f fl2|) for t (in particular solve: 
(t) = for t). Consequently we are unable to extract X Vl or Ai explicitly in terms of 
the parameters of the system. For brevity we omit the details and simply note that for 
the figures in this paper we determine X V1 and Ai by numerically solving transcendental 
expressions. This may be accomplished to any desired accuracy quickly and does not 
require the use of a differential equation solving method. 

The following two approximations are used in the analysis of later sections. For a 
wide range of parameter values the attracting periodic orbit passes close to the origin. If 
we approximate A„ 1 by finding where the next intersection of the forward orbit of (0,0) 
with the v-nullcline is (vi,wi), then we obtain 



A, 



071— SCT>7r ' 



(14) 



which is particularly accurate for e ~ e cr it, as shown in Fig. [3j 

Second, the two eigenvalues associated with {v* L ,w* L ) © are Pl iS i ow = 0(e) and 
PL.fast = Vl + 0(e). Within TZl, trajectories rapidly approach the associated slow eigen- 
vector. This eigenvector intersects the switching manifold, v — 0, at 



w L 



Xp 



L,slow 



a: 



(15) 



Consequently, trajectories such as large oscillations that spend a relatively long period 
of continuous time in TZl, exit this region extremely close to the point (0,wl)- This 
point is important below in the discussion of stochastic dynamics. It is usually suffi- 
cient to approximate Ai by considering (v^\t;0,w L ),w^(t;0,w L )) and the subsequent 



v 



(2), 



(t),w^ >(t)) and finding the value of A where (v ^ '(t), '(t)) intersects (1, 1). This 
is because Ai corresponds to the existence of a periodic orbit with a maximum v-value 
of 1, which must intersect (1, 1), see (j5]) and the surrounding discussion. 
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3 Exit distributions 



To analyze noise-driven MMOs we consider solutions to (pQ) with ([3]) in the presence of 
noise over long time frames such that transient behaviour has decayed. In this context 
we determine the fraction of oscillations that are small, the fraction that are medium, 
and the fraction that are large, referring to (jSj). One method is to simply solve the 
system for a long time and count the number of different oscillations. This Monte-Carlo 
approach is useful for obtaining a basic understanding of the system but poor for an 
accurate quantitative analysis because the system must be solved accurately for many 
parameter combinations requiring considerable computation time. Instead, since the 
system under consideration is PWL, we are able to use exit distributions for the regions 
(JTj) to approximate these fractions. This approach does not necessitate arbitrarily small 
e. In contrast, Muratov and Vanden-Eijnden [7] applied stochastic methods to (jTJ with 
([2]) by considering the system asymptotically (i.e. with arbitrarily small e and A) which 
essentially reduces the problem to one dimension. In [43J , the same system is considered 
but in the limit e — > which also reduces mathematical calculations to one dimension. 

Here we describe the exit distributions for forward orbits of ^ with (j3J) along various 
cross-sections of phase space. In the following section we use these exit distributions to 
identify MMOs. The four cross-sections we consider are: 



as depicted in Fig. [6j We exclude the switching manifold, v — 1, from calculations 
because large oscillations follow a sufficiently predictable path back to TZl when D <C 1. 

One method for computing a first exit distribution is to solve the Fokker-Planck 
equation for the probability density of the process (DQ) with (J3J) and absorbing boundary 
conditions (SI US] . Integration of the solution to this boundary value problem at the 
boundaries in an appropriate manner and over all positive time, may yield the desired 
exit distribution. However we dismiss this approach as it necessitates extensive numeri- 
cal computations, in part because drift dominates the diffusion which typically requires 
extra attention [H)J H7J US]. Instead we utilize the fact that within each region, IZj, the 
system is Ornstein-Uhlenbeck and, ignoring switching manifolds, has a known explicit 
solution |^5] . 

The transitional probability density, p[ , i.e. Pr((u(i), w(t)) G A | v(0) = vq,w(0) = 
w , (vq, Wo) £ Tli) =JJ Pt ( v i w \ v o, w o) dv dw, for the solution to the Ornstein-Uhlenbeck 



5 1 = {(0,w) 

5 2 = {{vi,w) 

£ 3 = {( v u w ) 
Z 4 = {{0,w) 



w < 0} U {(v, T]iV) I < v < vfi , 
W < W\} U {(v,T)iv) I V * < V < Vi} 
W > Wi} U {(VjTjiV) I V * < V < Vi} 

w > 0} U {(v, rjiv) I < v < v{} , 



(16) 



A 
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process of TZ±, i.e. ([I]) with f{y) = ijiv, after a time t is the Gaussian 



Pt(v,w \v ,w 



where 



27r v /det(6(t)) 



exp f-^A^e^A^J , 



(17) 



(18) 



w - u; (1) (t;t;o,Wo) 

The mean, (v^\ w^), is the solution to the system in absence of noise ( Fl2|) and 0(t) is 
the covariance matrix given by: 



G(t) = D 1 



Jo { e i2 ) ds f Q e 12 1 e 2 2 ds 
Jo e i2 e 22 ds J Q (e 2 2 ) ds 



D 



9u(t) e 12 {t) 
M*) e 22 {t) 



(19) 



where e^ lS denotes the (i, j)-component of the matrix exponential of A\s (jHJ) and we 
have introduced the 6ij for convenience. (Note that (1191) would contain more terms 
if ([T]) also included noise in the v equation.) The probability density (fT7|) obeys the 
Fokker-Planck equation 



dpt(v,w\v ,w ) 
dt 



where 



J t {v,w\v ,w ) 



-V ■ J t (v,«;|vo,tt;o) 
(771 u - 



e(m; — aw — X)p t 



d 2 apt 

2 dw 



(20) 
(21) 




v-nullcline 



exit 
distribution 



Figure 6: A sketch illustrating the exit distribution on E 2 (fIE|) for the forward evolution 
of a point on For clarity each Ej ( fl6|) is drawn with a different line type. 
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is the probability current [4"4"| |4"5] . By integrating ( 120]) and applying the divergence 
theorem, it follows that the net flow of probability across, say, the u-nullcline between 
(v*, u>l) and (v%, W\), is given by 



n ■ J t (v, rjiv\v ,w ) dv , 
where n is the normal vector of the f-nullcline pointing outwards [SJ US], i.e. here 

T 

. If trajectories were unable to cross the -u-nullcline more than once, 



n 



m 1 



then the integral 



n ■ J t (v, T)iv\v , w ) dt , (22) 







would be equal to the density of the first (and last) exit points for escape from below 
the u-nullcline. However, trajectories have multiple intersections with the f-nullcline 
due to the presence of noise. By considering two different time frames, we now show 
that these multiple intersections have a negligible effect and that (I2"2"j) represents an exit 
distribution suitable for our analysis. Specifically we first show that the probability of 
return to the f-nullcline after a short time is small. Then we show that within this 
short time frame points of multiple intersections are clustered. Finally we show that for 
longer time intervals after an intersection with the w-nullcline, trajectories are far from 
the nullcline, assuming small noise levels. 

We first look at return times for the i>-nullcline. Closed form expressions for first 
passage problems of mult i- dimensional Ornstein-Uhlenbeck processes are not straight- 
forward [HI [SO]; for this reason we simplify to a one- dimensional problem. Consider 
the forward orbit of a point on the i>-nullcline with v I < v < v%. Using y = w — rjxv to 
represent the distance from the nullcline, ([1]) with f(v) = rjxv may be written as 

dv = —y dt , 

/ \ (23) 

dy = {{Vi ~ a£ )y + (a — o'Vi)\ v ~ v i) £ ) dt + D dW , 



where we have substituted i>T = — - — . Intersections of the orbit with the nullcline are 

1 a— urn 



X 

-en 

determined by the y equation of (l2"3"j) which we conservatively reduce to 

dy = c dt + D dW , (24) 

where, for y > 0, the magnitude of the drift has a lower bound: 

c > (a- a?7i)(?j min - vl)e , 

assuming v > v min for some w min > For any 5 > 0, we are interested in Pr(y(t) = 
for some t > 5 \ y(0) = 0), i.e. the probability that a solution to ([24"]) with y(0) = 
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satisfies y(t) = at some t > S. To calculate this probability we let p(y, t) denote the 
transitional probability density for (|24|) with y(0) = 0, and condition over the event that 
y{8) = z, for all z e K: 

/oo 
p(z, 8) Pr(y(t) = for some t > 5 | y(S) = 
-oo 

Notice, 

Pr(y(t) = for some t > 5 | = 2) = Pr(y(t) = — z for some t > | ?/(0) = 0) , 
because (1241 has no explicit dependence on y and t. This enables us to write 



Pr(y(t) = for some t>6\ y(0) = 0) = / 

J — c 



p(z,5)G{-z)dz (25) 



where 

G(z) = Pr(y(t) = z, for some t > | y(0) = 0) . (26) 

G can be calculated from the density of the first hitting time of y(t) to z (refer to [51 
for more details) producing 



G{z) = c / p(z,t)dt . (27) 
Jo 

By using fl27|) and evaluating the integral on the right-hand side of fl25|) . we obtain 



Pr(y(t) = for some f > 5 \ y(0) = 0) = 1 - erf j • ( 28 ) 

For instance with D = 0.0012, e = 0.04 and (a, a) = (4, 1), whenever v min - v{ > 0.025 
the probability of return to the i>-nullcline after a time of 5 = 0.6 is less than 1%. 

Second, for the system (JTJ) with f{y) = r] X v we look at the distribution of future 
f-nullcline intersections up to a time 5. The solution ( TTTT) with wo = ijiVo evaluated on 
the w-nullcline and normalized is a Gaussian with mean and variance: 



-Ml2 ~ #22 + (^1^11 ~ 0l2 V W 

^22 - 2 Vl 9 12 + r^^n 



~ / x9 det 

U 2 2 — ll}\VYl. + ^i^ll 

respectively, where the 9ij were defined in ( TL91) . We observe that (1221) undergoes negligi- 
ble change when convolved by the Gaussian with ( 129|) and ( 130|) evaluated at t = 5. For 
this reason we use (122]) to compute exit distributions on the f-nullclines. The absence 
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of noise in the v equation of ([I]) ensures multiple rapid crossings through the switching 
manifolds are not permitted. Consequently we use an integral similar to (1221) for exit 
distributions across the other switching manifolds also. This accounts for all components 
of each £j (IT5|) . 

When the equilibrium, (v*,w*), is admissible, we expect the forward orbit of any 
point on Si to escape the lower half of TZ\ (below the v-nullcline) through £2. We 
calculate the exit distribution of the orbit through £2 with fl22|) . (Note, for simplicity 
we omit the term in J t ff2T]) when using f )22|) because it is dominated by the other 
terms in J t .) Using equally spaced data points and performing this calculation repeat- 
edly, we determine the exit distribution on £ 2 for any probability density of points on 
£1. From the exit distribution on £2 we continue in a similar fashion and compute 
the exit distributions on £3, £4 and lastly £1. Note these calculations use analytical 
expressions like (|22l) and not Monte-Carlo simulations. Numerically we observe that the 
iterative procedure of mapping a distribution on £1 to itself (through £2, £3 and £4) 
approaches the limiting distribution of the intersection of an arbitrary forward orbit of 
the system with £ 1; Fig. We use the limiting distributions on £ 2 and £ 3 to calculate 
the probability that an arbitrary oscillation is small, medium or large. The results for 
a range of parameter values are given in the next section. 

4 Mixed-mode oscillations 

In order to understand MMOs quantitatively, we say that ([I]) with ([3]) exhibits MMOs 
whenever both small and large oscillations occur at least 10% of the time. Specifically 
we find where exit densities corresponding to small and large oscillations both integrate 
to a value greater than 0.1. Fig. E] illustrates the dependence of MMOs on the primary 
bifurcation parameter, A, and the noise amplitude, D. Roughly the range of A values 
which permit MMOs increases with D. This matches our intuition, more noise allows for 
a wider variety of oscillations. We compute Fig. M using the iterative scheme described 
in £j3j Monte-Carlo simulations (not shown) give good agreement. 

MMOs exist in a region bounded on the left by the curve along which large oscillations 
occur 10% of the time and on the right by the curve along which small oscillations occur 
10% of the time. When D = the former curve has the value A = Ai, and the latter 
curve has the value A = X Vl . This is because the periodic orbit of the system in the 
absence of noise, ([H]), changes from small to medium at A = X Vl , and from medium to 
large at A = A a , §2J 

From Fig. |8l we see that MMOs do not occur for arbitrarily small D even near the 
canard explosion in contrast to what may be expected. This is because for D very 
small and X Vl < X < X\, medium oscillations dominate. Medium oscillations occur 
less frequently with increasing D. Note also that the MMO regions appear relatively 
symmetric with respect to A. 
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Figure 7: Stationary densities of intersections of (OQ) with (j3J) on Si (v — on the left, 
w = r/iv on the right). The solid curves in panel B are computed using the iterative 
method based on analytical expressions for the densities detailed in the text. The curves 
in panel A could also be calculated by this iterative prodecure, but instead it is more 
efficient to apply ( 1221) to the flow on the slow eigenvector of TZl with stationary variance. 
This is because in panel A rji is relatively large and oscillations enter TZl far from the 
origin and so are strongly attracted to the slow eigenvector of TZl- The histograms 
are calculated from a single trajectory of the system that was computed by numerical 
simulation over a time period of 2 x 10 5 . The value of wl ( |15l) is indicated in both 
panels. In panel A, w\ = 0.05 and A = 0.028; in panel B, w\ = 0.005 and A = 0.19. 
The remaining parameter values are D = 0.008, e = 0.04, V\ = 0.1, (a, a) = (4, 1) and 
(Vl,Vr) = (-2,-1). 



For the smooth system (CQ) with (T2]), we may roughly compute the region of MMOs, 
by the above definition, from Monte-Carlo simulations, Fig. |8]A. We see that MMOs 
occur over a smaller parameter range for the smooth version of the FHN model. This 
distinction is possibly explained by Fig. HI For the PWL model, all oscillations (including 
small oscillations) spend sufficient time in TZl to be drawn into the slow eigenvector of 
this region. Small and large oscillations are intertwined in TZl on their approach to TZ\] 
the amplitude of one oscillation is practically independent of the previous oscillation. 
(From a numerical viewpoint, in this situation fewer data points are required than in 
general.) In contrast, for the smooth system there is a significant distance between small 
and large oscillations and therefore more noise is required for, say, a large oscillation to 
follow a small oscillation. 

Noting this difference between the smooth and PWL models, we considered another 
parameter range for the PWL system that has a different exit distribution near the 
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Figure 8: Regions of MMOs defined by where at least 10% of oscillations are small and 
at least 10% are large. In panel A the parameter values used are the same as in Fig. |2j 
In panel B, W\ = 0.005; the remaining parameter values are unchanged. The region of 
MMOs for the smooth FHN model, ([T]) with (J2J), is superimposed in panel A. 



origin. For small values of the slope, rjx, still respecting ( FlUj) . the MMOs may include 
small oscillations that do not enter TZl, so that the exit distribution across Si may 
be bimodal, as shown in Fig. [Tj-B. Here large oscillations intersect Si near (0,wl) ( 1T51) 
whereas the majority of small oscillations intersect Si on the f-nullcline. We considered 
whether this type of bimodal exit distribution on Si plays a role analogous to distance 
between small and large oscillations in the smooth model, but we did not see any evidence 
of this effect. Specifically the MMO region, Fig. [8]-B, has a similar size and shape to 
the region in Fig. |8]A for which the corresponding value of r]i is an order of magnitude 
larger. 

The boundaries of the MMO regions shown in Fig. E] are relatively linear, hence we 
perform an analytical calculation of the slopes at D = 0. Let s sma ii [si arge ] denote the 
slope, at D = 0, of the curve along which 10% of oscillations are small [large]. 

Let us begin with the curve along which exactly 10% of oscillations are small. This 
curve intersects D = at A = X Vl at which the attracting periodic orbit created at A = 
intersects w = r/iv at v = v\. Here we can focus on small oscillations only, so it suffices 
to consider the linear systems oUZl and IZi, i.e. (JT]) with 



As A is increased, the maximum w-value of the deterministic periodic orbit increases at 




(31) 
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a rate, say, K\. Due to linearity, this rate is given simply by 

«i = -T- - (32) 

When A = A Vl , the periodic orbit intersects v = at some point (0, ) with 
tWA Bl < and the line w = r\\v at (vi,wi). If we now consider small D > but leave all 
other parameters unchanged, over a long time frame trajectories intersect v — at points 
approximately normally distributed about (0,10^). Since small oscillations neglect 
switching of (pQ) at v — v\, intersection points on w = rjiv are similarly approximately 
normally distributed about (vi,Wi), as shown in Fig. [HI with a standard deviation of 
say, 7i-D, where 71 is a constant that we compute below. The v-value of intersection 
points on w = r\\V then have the distribution N(vi + «i(A — A Vl ), 7i-D 2 ), using fl32|) . 
That is, if g sma ii denotes the probability density for these v-values, then 

n (,A - 1 -i^("-^i( A - A "i)) 2 ,ooN 

qsmalAV) - — == — e 7 i . [66) 



Then 10% of oscillations are small when 

%. m ^)*=^l-erf(-^l))=0,. (34) 

By rearranging the previous equation we deduce that the slope of the curve at D = is 

dD _ «4 

d\ V27! erf ^0.8) 

From §|2j K\ may be accurately calculated by solving transcendental equations. 

We obtain a good approximation to 71 as follows. Due to strong contraction in H L) 
the distribution of points on v = has a standard deviation that is much smaller than 
the standard deviation of points on w = r\\V. Than it is reasonable to approximate 
the distribution on v = by the single value w\ Vl • Then g sma n is equivalent to the exit 
distribution along the f-nullcline, thus by 



POO 

<? S maii(w) = w(v, T)iv) I pf> (v , 7]iv\0, W\ n ) dt + 0(D 2 ) , (36) 
Jo 



where w refers to (J6j). By ( TTTj) and (JHJ), 



ta («) = ™^/ ,„ } „ e-^dt + 0(D>), (37) 



2irD 2 J ^Q~Q, 



11^22 — C-12 



where 



2(^11^22 — & 



12/ 



M*) -e 12 (t) 
-M*) e n (t) 



v — 
rjiv — u/ 1 ) 



(38) 
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In the limit D — > 0, the asymptotic approximation to integral in (|37|) is determined 
from the main contribution of (p, which is its maximum value; formally this is achieved 
by Watson's lemma [53] • We omit the details of this calculation which produces 



7i = VKi • (39) 

We calculate the slope of the curve on which 10% of oscillations are large at D = in 
a similar fashion. Again we approximate the density of intersection points on v = by a 
point mass but this time we use the value wl ( 115)) and compute the density of intersection 
points on the switching manifold, v — V\. For small D this density is approximately 
Gaussian, i.e. N(w\ 1 + K2(X — Ai), ^\D 2 \ where when A = Ai the deterministic trajectory 
passes through the points (0, Wl) (or rather very near to this point), (v\, Wx^ and (1,1). 
If 9iar g e denotes this probability density, then 

-fA 1 -K2(A-Ai)) 2 



<7large(w) 



(40) 



'27T72-D 

The constant, k 2 , may be computed from f fl2|) using the chain rule for differentiation: 



K>2 



dt 



~8X + 



t=t in 



dw d\ dv* dX dwl dX 



where t int is the time taken for the trajectory to go from v 



dU 



dX 



+ 1 



+ 



dw dX dv\ dX dw\ dX 




to v 
d v (i) 



V\ and 



dt 



(41) 



(42) 



t=t in 




0.04 0.08 v. 0.12 0.16 



Figure 9: Intersections of ([T]) with fl3Tl) 01110 = rjiv using the same parameter values as 
Fig. [2] with also A = 0.028, D = 0.0004. The probability density curve is given by ( |33|) 
and the histogram is calculated from a single numerically computed trajectory solved 
up to a time 2 x 10 5 (part of which is shown also). 



18 



By a calculation similar to that for 71 described above, we obtain 




72 = W 011 (j) - 2012 (j) + 022 , (43) 

using (jSJ). 

Unlike small oscillations, large oscillations traverse TZ 2 and so we must also 
consider the flow in these regions. For D = and A near Aj, we let ^2 denote the first 
intersection of the backwards orbit from (1,1) with v — v± so that we may distinguish 
medium and large oscillations on v — V\ when D = 0. We write 

w 2 = w Xl + k s (X - Xi) , (44) 

ignoring higher order terms and where % may be calculated in a manner similar to k 2 . 
For small D > 0, the forward orbit of any point (v\, w), with w < w%, has the probability, 
Pia.rge(w, A, D) , of undergoing a large, rather than medium, oscillation before returning 
to IZl- We find that for small D there is a very sharp transition of p\ aTge at w = w 2 so 
that it suffices to use the approximation pi arge (u>) = H{w 2 — w), where H(z) = for 
z < and H(z) = 1 for z > 0. Consequently, 10% of oscillations are large when 

: + ^^-)---K i --( fiiz ^))-"- <45) 

where gi arg e is given by pUj) . By rearranging this expression we arrive at 

dD _ k 2 - K3 f s 

Slarge " dX ~ V2 l2 erf-^O.S) ' 1 J 



We have verified that computation of s sma u and si arge by the expressions ( 1351) and ( 1461 
matches with computation by exit distributions ( 1221) as in Fig. [HJ 

The variation of the slopes ( 1351) and ( 146|) with respect to 771, e and a is shown in 
Fig. ED Both (|35|) and ( 146|) approach zero as e — > 0, but accounting for the fact that the 
noise amplitude, D, is not multiplied by e in (JT]), the scaled values ^s sma u and ^si arge vary 
relatively slightly. With increasing iji , s ama ii an d — -Siarge increase slightly and approach 
the same value; with increasing a, s sma n and — si arge decrease slightly. Due to linearity, if 
T]i is held constant and V\ is increased, s S maii is unchanged. We do not need to consider 
variation in wi and 772 because these values may be written in terms of rji and V\ (j4|). 
Additionally, the slopes are not strongly affected by i] L (as long as r/ L 0) and i] R . 
Therefore, for the most part, |s sma u and — ^si arge lie in, say, the interval [1.8, 2.5]. Hence 
for intermediate values of although the interval of A values which permit MMOs varies 
widely with the system parameters, the width of this interval is robust with respect to 
parameter change. For small values of — , MMOs may not occur at all. For very large 



values of — , multiple crossings on the Sj may generate different results. 
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Figure 10: The dependence of X V1 , X\, s S maii and si arge on parameter values. The lower 
curves are X Vl and Ai; the upper curves are the slopes. The X Vl curves have end points 
at (771, A) = (e<7, ~(a — ea 2 )vi) and (2y/ea — ea, 0). In panel A, a = 4. In panel B, 
e = 0.04. In both panels, a — 1, t>i = 0.1, and (771,, 77^) = (—2, —1). We have scaled the 
slopes by ~ because the noise amplitude D is not multiplied by e in (CQ). 

5 Conclusions 

We have studied MMOs in a PWL version of the FHN model, (DQ) with (EJ). To ob- 
tain quantitative results we have defined oscillations as small, medium, or large by the 
maximum v- value attained §5§. Furthermore we define MMOs by where at least 10% of 
oscillations are small and at least 10% are large (this approach may be applied to any 
preferred values of the percentages). We incorporate noise additively in one equation 
for transparency of analysis. Numerically we have observed that additive noise in both 
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equations yields similar mixed-mode dynamics. 

The existence of a canard explosion in the system with no noise is a consequence 
of incorporating a four-piece PWL function into the model. We have introduced an 
analogy for canard points of smooth systems, specifically we identify two values, X Vl 
and Ai, at which the periodic orbit of the deterministic system changes from small to 
medium, and from medium to large. Both values increase with increasing a, and e, and 
decrease with increasing rji as shown in Fig. [TUJ 

Near the canard explosion noise drives MMOs. The boundaries of the regions of 
MMOs shown in Fig. |H] are approximately linear for small D and we have calculated 
their slopes, at D — 0. For small values of D we have shown that the probability 
current provides an approximation for exit distributions; for large values of D this 
approximation may no longer be valid. Unless the noise amplitude is extremely small, 
MMOs exist over some interval of A values. We have illustrated that for constant — 

e 

typically the width of this interval changes minimally with a relatively large variation 
in the values of the other system parameters. 

For the results in this paper we have used the value t]l = —2 for the slope of the 
f-nullcline for v < 0. With similar or more negative values of t)l the system exhibits the 
same qualitative behaviour. However, with larger values of r/ L , say — 1 < 77^ < 0, ([TJ) 
with ([3]) may have multiple attracting solutions in the absence of noise. The coexistence 
of attracting small and large periodic orbits, naturally produces MMOs in the presence 
of noise though this is via a different mechanism than the one studied here. 

The FHN model of Makarov et. al. [8J, which includes additional nonlinearity, read- 
ily exhibits MMOs. It is possible that this addition strengthens the attraction of the 
periodic orbit for v < 0, mimicking the slow eigenvector discussed above and causing 
MMOs to be more robust. 
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